
#!usr/bin/perl -w
use strict;

#this script is to extract RT domains from ltrharvest result fasta file
# perl extract_blastx.pl -db Ncra.ltrharvest.fa -p Ncra.ltrharvest.fa.blastx.tab
use Bio::DB::Fasta;
use Bio::SeqIO;
use Getopt::Long;

my $db;
my $prefix = './domains/';
GetOptions(
           'd|db:s' => \$db,
           'p|prefix:s' => \$prefix,
           );

mkdir($prefix) unless -d $prefix;
if( ! defined $db ) {
    die"Must provide a LTR db with -d or --db option\n";
}

my $dbh = Bio::DB::Fasta->new($db);
my %fams;

while(<>) {
#	my @word=split(/\s+/,$_);	
#	my ($seqid,$start,$end) = ($word[0], $word[6],$word[7]); # get seq name, start and end from blastx table

    # another way to do this is to assign all the data that comes from the list
    # if you can't remember the order, run a test blast (1 seq against a db)
    # and use the -m 9 option (instead of -m 8). It will print a header
    # listing the names of the columns
    my ($query,$target,$percent_id, $aln_length, $mismatches, $gaps,
	$qstart,$qend, $hstart,$hend, $evalue, $bitscore) = split(/\s+/,$_);
    if( ! exists $fams{$query} ) {
	# do you want to only get the best hit for the query?
	# just store one value and don't do the 
	# 
	# push @{$fams{$query}}, [$qstart,$qend,$target,$hstart,$hend];
	$fams{$query} = [$qstart,$qend,$target,$hstart,$hend];
	
    } else {
	# we are skipping all the rest of the hits for this query
	# - this is based on E-Value since that is how the hits are sorted, but
	# we may want to invoke another criteria like length of alignment?
	# that would require keeping the longest alignment 	
    }
}

my $out = Bio::SeqIO->new(-format => 'fasta', 
			  -file   => ">$prefix/blastx_domains.fas");

for my $ltrname (keys %fams) {    
    
    my ($qstart,$qend,$target,$hstart,$hend) = @{$fams{$ltrname}};
    
    my $seq = $dbh->seq($ltrname,$qstart=>$qend);
    my $seqobject = Bio::PrimarySeq->new(
					 -id   => "$ltrname\_$qstart\_$qend",
					 -desc => "hit to $target:$hstart..$hend", 
					 -seq  => $seq);
    $out->write_seq($seqobject);
}

